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Abstract 

We report single-cluster Monte Carlo simulations of the Ising model 
on three-dimensional Poissonian random lattices with up to 128 000 ~ 
50 3 sites which are linked together according to the Voronoi/Delaunay 
prescription. For each lattice size quenched averages are performed 
over 96 realizations. By using reweighting techniques and finite-size 
scaling analyses we investigate the critical properties of the model in 
the close vicinity of the phase transition point. Our random lattice 
data provide strong evidence that, for the available system sizes, the 
| resulting effective critical exponents are indistinguishable from recent 

high-precision estimates obtained in Monte Carlo studies of the Ising 
model and </> 4 field theory on three-dimensional regular cubic lattices. 



1 Introduction 



Experimental studies of the critical behavior of real materials are often con- 
fronted with the influence of impurities and inhomogeneities. For a proper 
interpretation of the measurements it is, therefore, important to develop a 
firm theoretical understanding of the effect of such random perturbations. In 
many situations the typical time scale of thermal fluctuations in the ideal- 
ized "pure" system is clearly separated from the time scale of the impurity 
dynamics, such that to a very good approximation the impurities can be 
treated as quenched (frozen), random disorder. 

The importance of the effect of quenched, random disorder on the critical 
behavior of a physical system can be quite generally classified by the critical 
exponent of the specific heat of the pure system, a p . The Harris criterion 1 
asserts that for a p > quenched, random disorder is a relevant perturbation, 
leading to a different critical behavior than in the pure case. In particular one 
expects 2 in the disordered system that v > 2/D, where v is the correlation 
length exponent and D the dimension of the system. Assuming hyperscaling 
to be valid, this implies a = 2 — Du < 0. For a p < disorder is irrelevant, 
and in the marginal case a p = no prediction can be made. For the case of 
(non-critical) first-order phase transitions it is known that the influence of 
quenched, random disorder can lead to a softening of the transition. 3 

Many theoretical and numerical studies have been devoted to quenched, 
random site-dilution (SD), bond-dilution (BD), or more general random-bond 
(RB) systems. Since for the three-dimensional (3D) Ising model it is well 
known that a p ~ 0.1 > 0, quenched, random disorder should be relevant for 
this model. This has indeed been verified by a variety of techniques: Monte 
Carlo (MC) simulations for SD 4 ~ 6 and BD, 7 ' 8 high-temperature series (HTS) 
expansions for BD 9 and field theoretical renormalization group studies. 10-13 
For an excellent review and an extensive list of experimental, theoretical and 
numerical estimates in the last two decades, see Ref.. 14 As a result consensus 
has been reached that, while the critical exponent ratios 7/V = 2 — 77 and 
Pjv are hardly distinguishable from the pure model, the correlation length 
exponent v clearly signals the disordered fixed point, 

7/1/ = 1.966(6), (3/v = 0.517(3), v = 0.6304(13) (pure 15 ), 

7/1/ = 1.970(3), p/v = 0.515(2), v = 0.678(10) (disordered 13 ), 

and, moreover, satisfies the bound z/>2/.D = 2/3in the disordered case. 
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Recently also the predicted softening effect at first-order phase transitions 
has been confirmed for 3D g-state Potts models with q > 3 using MC 16-18 
and HTS 19 techniques. The overall picture is even better in two dimensions 
(2D) where several models with a p > (SD Baxter model, 20 SD Baxter- 
Wu model, 21 3- and 4-state RB Potts model, 22 Ashkin- Teller model 23 ) and 
the marginal (a p = 0) 2D Ising model 24-28 have been investigated. Also the 
softening of the first-order transitions for the 2D models with q > 5 has been 
confirmed. 29-32 For a recent review, see Ref.. 33 

In this paper we study another type of quenched, random disorder, namely 
connectivity disorder, a generic property of random lattices whose local coor- 
dination numbers vary randomly from site to site. Physically the concept of 
random lattices plays an important role in an idealized description of the sta- 
tistical geometry of random packings of particles. 34-36 A prominent example 
is the crystallization process in liquids, and many statistical properties of ran- 
dom lattices have been studied in this context. 37 From a more technical point 
of view, random lattices provide a convenient tool to discretize spaces of non- 
trivial topology without introducing defects or any kind of anisotropy. 38-40 
This desirable property has been exploited in a great variety of fields, ranging 
from diffusion limited aggregation 41 and growth models for sandpiles 42 over 
the statistical mechanics of membranes and strings 43 to quantum field theory 
and quantum gravity. 38-40 ' 44 ' 45 The preserved rotational, or more generally 
Poincare, invariance suggests that spin systems or field theories defined on 
random lattices should reach the infinite-volume or continuum limit faster 
than on regular lattices. Conceptually, however, such an approach does only 
make sense as long as the critical properties of the considered system are 
not modified by the irregular lattice structure. In view of the quite general 
Harris criterion this is a non-trivial question (in particular due to the in- 
herent spatial correlations of the disorder in this case), at least for systems 
characterized by a p > 0. 

Specifically we considered 3D Poissonian random lattices of Voronoi/De- 
launay type, and performed an extensive computer simulation study of the 
Ising model for lattices varying in size from N = 2 000 ~ 13 3 to 128 000 ~ 50 3 
sites. For each system size quenched averages over the connectivity disorder 
are approximated by averaging over 96 independent realizations. We con- 
centrated on the close vicinity of the transition point and applied finite-size 
scaling (FSS) techniques 46 to extract the critical exponents and the (weakly 
universal) "renormalized charges" £/| and U\. To achieve the desired accu- 
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racy of the data in reasonable computer time we applied the single-cluster 
algorithm 47 ' 48 to update the spins and furthermore made extensively use of 
the reweighting technique. 49 

Previous studies of connectivity disorder focused mainly on 2D where pro- 
nounced effects were observed in MC simulations of g-state Potts models on 
quenched, random graphs provided by models of quantum gravity (modified 
universality classes for q = 2 and 4, 50-52 and softening for q = 10 52 ' 53 ). For 
2D random lattices of Voronoi/Delaunay type, on the other hand, no influ- 
ence was seen in simulations for q = 2 54 ' 55 and q = 8. 56 The main difference 
between these two types of random lattices is the highly fractal structure 57 
of random gravity graphs which suggests a stronger "degree of disorder". 
A similar dependence on the "degree of disorder" was recently observed for 
several aperiodic perturbations. 58 

The rest of the paper is organized as follows. In Sect. 2, we describe some 
properties of the random lattices used in the simulations, and in Sect. 3 
we define the model and give a few simulation details including estimates 
of autocorrelation times. The quantities measured are defined in Sect. 4, 
where also their theoretically expected FSS behavior is recalled. In Sect. 5, 
we present the FSS analysis of our data close to the transition point which 
yields estimates for the critical exponents. Finally, in Sect. 6 we present a 
discussion of our main results and close with a few concluding remarks. 



2 Random lattices 

A. Voronoi/Delaunay random lattices: The notion of "a random lat- 
tice" is by no means unique and, in fact, many different types of random 
lattices have been considered in the recent literature. 41 ' 50 ' 52 ' 53 ' 59 In this pa- 
per we concentrate on so-called Poissonian Voronoi/Delaunay random lattices 
which, in arbitrary dimensions, are defined as follows 38 ' 44 ' 45 : 

1. Draw N sites Xi at random in a unit volume (square in 2D, cube in 3D, 
...). 

2. Associate with each site X{ a Voronoi cell, q = {x \ d(x,Xi) < d(x,Xj) 
Vj 7^ i}, which consists of all points x that are closer to the center 
site Xi than to any other site. Here d(x, y) denotes the usual Euclidean 
distance. This yields an irregular tessellation of the unit volume with 
D-dimensional Voronoi cells (polygons in 2D, polyhedra in 3D, . . . ). 
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3. Construct the dual Delaunay lattice by linking the center sites of all 
Voronoi cells which share a common face. 

The first step approximates the Poissonian process. In the second step we al- 
ways assume periodic boundary conditions, i.e., construct lattices of toroidal 
topology. Using this construction the number of nearest neighbors of the 
Delaunay lattice, the local coordination number q, varies randomly from site 
to site. This constitutes the special type of quenched, random connectivity 
disorder we are investigating in this work. 

The actually employed construction of the random lattices follows loosely 
the method of Ref. 60 and is described in the Appendix. Following this method 
we succeeded to reduce the complexity of the lattice construction for all 
practical purposes from order N 2 , as expected for the most straightforward 
implementation, to order N, up to a small overhead oc iV 2 resulting from the 
initial calculation of the distances between all pairs of two sites. The actually 
measured CPU times as a function of iV are shown in Fig. I. 

B. Random lattice properties: To test our random lattice construction 
we measured several quantities which characterize the topology and geometry 
of the lattices. These quantities are exactly known in the limit N — > oo. 
"Topological" quantities are the number of links iVj, the number of triangles 
and the number of tetrahedra N T , which according to Euler's theorem 
should satisfy for a torus topology in any number of dimensions and for any 
number of sites N the relation 

N-N l + N A -N T = 0. (1) 

We have of course explicitly tested that this topological constraint is satisfied 
for all realizations, which is quite a sensitive confirmation that the lattice 
construction works properly. The measured averages of q, Ni, and N T 
over the 96 realizations used in the simulations are collected in Table 1. The 
error bars are computed from the fluctuations among the 96 realizations. We 
see that the analytically known N — > oo limits are approached very rapidly. 
The only exception is perhaps our smallest lattice with N = 2 000 sites where 
deviations from the infinite- volume limit are clearly visible. For N > 4 000 
in particular the average number of nearest neighbors, q, is fully consistent 
with the theoretical value of 

48 9 

g = 2 + —7r 2 = 15.5354.... (2) 
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Table 1: The average coordination number q and the total number of simplices 
normalized by the number of sites. The error bars are computed from the 
fluctuations among the 96 realizations. 



N 

2 000 
4 000 
8 000 
16 000 
32 000 
64000 
128 000 
exact 
exact 



Q 

15.544(4) 

15.532(3) 

15.537(2) 

15.535(2) 

15.534(1) 

15.5351(7) 

15.5349(5) 

15.5354. . . 

2 + I- 2 



7.7721(19) 
7.7661(16) 
7.7683(10) 
7.76744(76) 
7.76706(55) 
7.76755(36) 
7.76743(26) 
7.76772. . . 
1 + "' 



35 



7T- 



N A /N 

13.5441(37) 

13.5321(32) 

13.5366(21) 

13.5348(15) 

13.5342(11) 

13.53507(72) 

13.53486(53) 

13.53545. . . 

35" 



N T /N 

6.7721(19) 

6.7661(16) 

6.7683(10) 

6.76744(76) 

6.76706(55) 

6.76755(36) 

6.76743(26) 

6.76772. . . 

35 " 



As "geometric" quantities we measured the average volumes of the sim- 
plices, i.e., the average link length I = Y^Zi ((J2ie T h)/6) /N T , the average 
triangle area A = (£•=! Ai)/N A , and the average volume of a tetrahedron 
t — (J2i^iT~i)/N T . Notice that the average link length is defined in such a 
way that we first average over the sides of a given tetrahedron and then over 
all the tetrahedra. This is thus the average link length per tetrahedron, and 
not the mean link length averaged over the whole lattice, Z (A ° = (£^i k)/N h 
which is always larger due to the fluctuations in the number of tetrahedra 
from realization to realization. Our results normalized with an appropriate 
power of the density p = N/V are displayed in Table 2. From N = 4 000 
on the numerical results are again fully consistent with the analytical pre- 
dictions as iV — > oo. Since all numbers agree very well with the analytical 
predictions we can be quite sure that our lattice construction works properly 
and that we have picked a typical sample of random lattices. 

The (normalized) probability densities P(q) of the coordination numbers 
q are shown in Fig. 2 for the lattices with N = 64 000 and N = 128 000 

sites. The average coordination number is q = 2 + ||7r 2 = 15.5354 Due 

to the long tail of P(q) for large values of q, the maximum of P(q) is found 
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Table 2: Average simplex volumes properly normalized to natural units. 



N 



7 (A V 1/3 



l/p 



-1/3 



A/p 



-2/3 



r/p- 



2 000 
4000 
8000 
16 000 
32 000 
64 000 
128 000 
exact 
exact 



1.28540(22) 

1.28566(18) 

1.28569(13) 

1.285623(75) 

1.285554(68) 

1.285537(40) 

1.285541(28) 



1.23670(17) 

1.23717(13) 

1.23711(10) 

1.237091(55) 

1.237027(49) 

1.237002(32) 

1.237012(20) 

1.237033. . . 

l 4vr J 3 2 -4 4i ^3^ 



0.597004(95) 
0.597362(70) 
0.597312(55) 
0.597306(31) 
0.597289(24) 
0.597272(19) 
0.597275(11) 
0.597286. . . 

\i-K> 3 s 7T 1 V3/ 



0.147666(41) 

0.147797(34) 

0.147747(23) 

0.147766(17) 

0.147775(12) 

0.1477630(78) 

0.1477666(58) 

0.1477600... 

35 
24tt 2 



at the next smaller integer number q 
of 12.03%. 



15, which occurs with a probability 



3 The model and simulation techniques 

A. Model: We simulated the standard partition function of the Ising model, 
Z = Te- KE ; E = -Ts lSi ; Si = ±l, (3) 



{*} 



^ ' SjSj, Si 
(ij) 



where K = J/ksT > is the inverse temperature in natural units, and (ij) 
denotes the nearest-neighbor links of our three-dimensional toroidal random 
lattices. In (3) we have adopted the convention used in Ref. 55 and assigned 
to each link the same weight. 

Another interesting option would be to assign to each link a weight de- 
pending on the geometrical properties of the Voronoi/Delaunay construction 
such as the length of the links or suitably defined areas of the associated tes- 
selation. In addition to the connectivity disorder this would introduce also a 
(correlated) disorder of random-bond type. In order not to mix up these two 
quenched disorder types, we decided to concentrate in this study exclusively 
on the effect of the connectivity disorder. 
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B. Simulation: The finite-size scaling (FSS) analysis is performed on the 
basis of seven different lattice sizes with N = 2 000, 4 000, 8 000, 16 000, 
32 000, 64 000, and 128 000 sites. For later use we adopt the notation for 
regular lattices and define a linear lattice size L by L = iV 1 / 3 . The linear 
sizes of the lattices thus vary from L fa 12.6 to L fa 50.4. To investigate the 
dependence of thermal averages on different realizations we performed MC 
simulations, for each lattice size, on the 96 randomly chosen random lattice 
realizations discussed in Sect. 2. 

For the update of the Ising spins we employed Wolff's single-cluster (1C) 
algorithm. 47 Various tests, in particular for the Ising model on two- and 
three-dimensional regular lattices, have clearly demonstrated that critical 
slowing down can be significantly reduced with this non-local update al- 
gorithm. 61-63 These tests also showed that in particular in three dimensions 
the single-cluster algorithm is more efficient than the original multiple-cluster 
formulation of Swendsen and Wang. 48 

All runs were performed in the vicinity of the critical point K c . As a 
first rough guess of K c we used the mean- field bound K c > Kmf = l/<? ~ 
0.064. Due to the large average coordination number of 3D random lattices 
the mean-field approximation should be better than for the simple cubic 
(SC) lattice where Kfp' = 1/6 is about 1.33 times smaller than the actual 
}{(sc) ^ q 221 654 6. Therefore, by applying the same correction factor to the 
random lattice mean-field estimate we expect that K c is bounded from above 
by K c < 0.085. This heuristic argument thus suggests that 0.064 < K c < 
0.085, such that K c fa 0.075 should be a reasonable a priori guess. Once 
good simulation points Kq were estimated for the two smallest lattices with 
N = 2 000 and 4 000 sites by determining for a few realizations the location 
of the specific-heat and susceptibility maxima, we used FSS extrapolations 
(assuming v fa 0.63) to predict K for the larger lattices. 

Depending on the lattice size from 30 000 to 180 000 clusters were dis- 
carded to reach equilibrium from an initially completely ordered state. Pri- 
mary observables are the energy per spin, e = E/N, and the magnetization 
per spin, m = M/N = J2i$i/N, which were measured every nmp cluster 
flip and recorded in a time-series file. The average cluster size (\C\) is 
an estimator for the reduced susceptibility in the high-temperature phase, 
Xred = N(m 2 ), and therefore scales with the lattice size according to L 7 ^, 
where 7 and v are the usual critical exponents. Since in three dimensions 
7/V = 2 — i] fa 2, we have to perform oc L cluster flips in order to flip, on 
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Table 3: Monte Carlo parameters of the simulations. N = L 3 is the lat- 
tice size, the hash mark symbol # denotes the number of realizations, K 
the inverse simulation temperature, n t hcrm the number of cluster flips during 
equilibration, and n mcas the number of measurements taken every 7ifl ip cluster 
flip. 



N 


L 


# 


K 


"therm 


^meas 


"-flip 


2 000 


12.6 


96 


0.0735 


30 000 


100 000 


4 


4000 


15.9 


96 


0.0735 


50 000 


100 000 


5 


8 000 


20.0 


96 


0.0732 


66 000 


110000 


6 


16 000 


25.2 


96 


0.0729 


80 000 


100 000 


8 


32 000 


31.7 


96 


0.0728 5 


100 000 


100 000 


10 


64000 


39.1 


96 


0.0726 5 


150 000 


100 000 


13 


128 000 


50.4 


32 


0.0725 3 


180 000 


100 000 


16 


128 000 


50.4 


64 


0.0725 9 


180 000 


100 000 


16 



the average, approximately the same fraction of the total number of spins, 
N = L 3 , for all lattice sizes. By adjusting the absolute scale of we made 
sure that for all lattice sizes the measurements were taken after about N/2 
spin flips. Since for the single-cluster update algorithm the autocorrelation 
times scale only weakly with L (see the discussion below) one then expects 
that, with about the same number of measurements, the statistical accuracy 
is comparable for all lattice sizes. With this set-up we performed 100 000 
measurements for each lattice size and realization (110 000 for iV = 8 000). 
For more details on the employed statistics, see Table 3. 

C. Update dynamics: A useful measure of the update dynamics is the 
integrated autocorrelation time f. 64 To estimate f for the measurements 
of the energy e and the magnetization m we used the fact that f enters in 
the error estimate e 2 = o" 2 2f/n meas for the mean value O of n mcas correlated 
measurements with variance a 2 = (O; O) = (O 2 ) — (O) 2 , and determined e 2 
by blocking procedures. Using 100 blocks of 1000 measurements each we ob- 
tained the results compiled in Table 4 where all results are averaged over the 
96 randomly chosen realizations. We see that the integrated autocorrelation 
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Table 4: Average cluster size [(|C|)] av and autocorrelation times of energy 
and magnetization at the simulation point K , where t = f ■ f and f = 
7ifli p (|C|)/iV. The tau's are obtained from a blocking analysis on the basis of 
100 block. 



N 


[<|C|}]av 


[^e ] av 




[/lav ■ [r e 


av [^m ] av 


[^~m] av 


[/]av ' [^m]av 


2 000 


246.9(1.9) 


2.529(40) 


1.245(19) 


1.249 


2.434(38) 


1.200(20) 


1.202 


4000 


453.3(3.9) 


2.443(42) 


1.375(20) 


1.384 


2.412(38) 


1.360(20) 


1.367 


8 000 


740.1(5.5) 


2.702(42) 


1.494(22) 


1.500 


2.594(40) 


1.436(22) 


1.440 


16 000 


1084.0(9.5) 


2.955(43) 


1.593(21) 


1.602 


2.730(38) 


1.475(22) 


1.480 


32 000 


1962(17) 


2.741(50) 


1.666(24) 


1.681 


2.572(44) 


1.567(23) 


1.577 


64 000 


2678(24) 


3.177(55) 


1.718(27) 


1.728 


2.795(42) 


1.516(23) 


1.520 


128 000 


3521(87) 


3.92(11) 


1.705(42) 


1.725 


3.157(72) 


1.385(42) 


1.390 


128 000 


4325(41) 


3.399(72) 


1.824(31) 


1.838 


2.964(60) 


1.593(29) 


1.602 



times for the measurements of e and m are of the order of f e ~ 2.5 — 3.5 and 
T m 2.4 — 3.0. Since completely uncorrelated data correspond to f = 0.5, 
our thermal sample thus effectively consists of about 15 000 — 20 000 uncorre- 
lated measurements for each of the 96 realizations. This amounts to a total 
of (1.5 — 2.0) x 10 6 effectively uncorrelated measurements for each lattice 
size. 

While this properly characterizes the effective statistics of our simula- 
tions, the numbers for f of a single-cluster simulation are not yet well suited 
for a comparison with other update algorithms or with single-cluster simu- 
lations on regular lattices. To get a comparative work-estimate, the usual 
procedure 47 is to convert the f by multiplying with a factor / = nm. P {\C\)/N 
to a scale where, on the average, measurements are taken after every spin has 
been flipped once (similar to, e.g., Metropolis simulations). For quenched, 
random systems this procedure is not unique due to the necessary average 
over realizations. In Table 4 we therefore present both options, [r] av = [/-T] av 
and also [/] av • [f] av . The differences between the two averaging prescriptions 
are, however, extremely small. 

The numbers in Table 4 obtained in this way are very similar to results 
for the regular simple cubic (SC) lattice. 61,62 By fitting a power law, r oc L z , 
to the data for the five largest lattices we obtain [r e ] av = 0.82(6)L a20 ( 2 ), and 
[Tm]av = 1.07(8)L ai0 ( 2 \ respectively. The dynamic critical exponents z for 
the random lattice simulations should be compared with z e = 0.28(2) and 
z x = 0.14(2) for the SC lattice. 61 
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4 Observables and finite-size scaling 



From the time series of e and m it is straightforward to compute in the FSS 
region various quantities at nearby values of Kq by standard reweighting 
methods. 49 For the estimation of the statistical (thermal) errors, for each of 
the 96 realizations the time-series data was split into 100 bins, which were 
jack-knifed 65 to decrease the bias in the analysis of reweighted data. The final 
values are averages over the 96 realizations which will be denoted by square 
brackets [. . .] av , and the error bars are computed from the fluctuations among 
the realizations. Note that these errors contain both, the average thermal 
error for a given realization and the theoretical variance for infinitely accurate 
thermal averages which is caused by the variation of the quenched, random 
geometry of the 96 lattices. 

From the time series of the energy measurements we can compute by 
reweighting the average energy, specific heat, and energetic fourth-order pa- 
rameter, 



u(K) 
C(K) 

V(K) 



l(E)UN, 
iT 2 iV[(e 2 )-(e) 2 ] 

[1- ^ 



3(e 2 ) 2J 



(4) 
(5) 

(6) 



Similarly we can derive from the magnetization measurements the average 
magnetization, the susceptibility, and the magnetic cumulants, 



m{K) = [(H)]av, 

X (K) = K N[(m 2 ) - 

U 2 (K) = [1- <™ 2 > 

U 4 (K) = [1- 



\m\y 



3(H) 2jav ' 
(m 4 ) 



(7) 
(8) 

(9) 
(10) 



3(m 2 ) 2 

Further useful quantities involving both the energy and magnetization are 
d[(\m\)U 



dK 
dln[(\m\) 
dK 



\m\E)-(\m\)(E)l 
(HE) 



(H) 



(E) 



(11) 
(12) 
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d\n[(m% 
dK 



(m 2 E) 
(m 2 ) 



-(E) 



(13) 



In the infinite-volume limit these quantities exhibit singularities at the 
transition point. In finite systems the singularities are smeared out and scale 
in the critical region according to 



C 

[<H>W 

X 

d[(M)U 
dK 
d\n[(\m\P)] ay 



dK 



dU, 



2p 



dK 



C rcg + L^f c (x)[l + .. 
L-^f m (x)[l + ...], 
I?f"f x (x)[l + ...], 

L^^f m/ (x)[l + ...], 

lVv p (x)[i + ...], 

L^KX*) [! + •••], 



(14) 
(15) 
(16) 

(17) 
(18) 
(19) 



where C rcg is a regular background term, u, a, (3, and 7, are the usual critical 
exponents, fi(x) are FSS functions with 



x — (K — K c )L l/u 



(20) 



being the scaling variable, and the brackets [1 + . . .] indicate corrections-to- 
scaling terms which become unimportant for sufficiently large system sizes 
L. 



5 Results 

By applying standard reweighting techniques to each of the 96 time-series 
data we first determined the temperature dependence of Ci(K), Xi{K), . . . , 
% — 1, . . . , 96, in the neighborhood of the simulation point Kq. By estimating 
the valid reweighting range for each of the realizations we made sure that 
no systematic errors crept in by this procedure. Once the temperature de- 
pendence is known for each realization, we can easily compute the disorder 
average, e.g., C(K) = Ei=i Q (K) /96, and then determine the maxima of 
the averaged quantities, e.g., C max (-ft' maxc ) = max^ C(K). The locations of 
the maxima of C, x, dU 2 /dK, dU^/dK, d[(|m])] av /d£T, dln[(|m|)] av /(Z.K', and 
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dln[{m 2 )] av /dK provide us with seven sequences of pseudo-transition points 
if mEHi (L) which all should scale according to K maXi (L) = K c + diL^ 1 ^ + . . .. 
In other words, the scaling variable x = (K max .(L) — K C )L 1 / U = Oj + . . . 
should be constant, if we neglect the small higher-order corrections indicated 
by . . .. To give an idea on how these sequences approach K c we show below 
in Table 6 besides the estimates for K c also the values of Oj, i — I, ... ,7 (see 
also Fig. 4). 

It should be emphasized that while the precise estimates of Oj do depend 
on the value of u, the qualitative conclusion that x ~ const, for K maXi does 
not require any a priori knowledge of v or K c . Using this information we 
have thus several possibilities to extract unbiased estimates of the critical 
exponents u, a/u, f3/v, and 7/1/ from least-squares fits assuming the FSS 
behaviors (14)-(19). Once v is estimated we can then use K max .(L) = K c + 
OiL -1 /" + ... to extract also K c and Oj. 

A. Critical exponent v. Let us thus start with the correlation length 
exponent v for which we can obtain 4 x 7 = 28 different estimates by consid- 
ering the scaling of drn[(|m|)] av /d£T, <iln[(m 2 )] av / dK , dU 2 /dK, and dU^/dK 
at the seven sequences of pseudo-transition points K maXi (L). Of course, these 
estimates are statistically not uncorrelated, but they are differently affected 
by corrections to the leading FSS behavior. To test the importance of these 
corrections to scaling we first estimated v from fits using all available lattice 
sizes, and then compared with the results of fits where the two smallest sizes 
were successively discarded. As a result we did not observe any systematic 
improvement by omitting the smallest lattices. In fact, already for the fits 
using all sizes, we obtained goodness-of-fit parameters 66 Q > 0.3 for about 
80% of the 28 fits. The only unacceptable fit was that of d\n[(m 2 )] av /dK 
at the K mSLXx sequence with Q = 0.01. Here we omitted the iV = 2000 data 
point which improves the goodness to Q — 0.22. 

This analysis clearly shows that with our data there is no need to include 
corrections-to-scaling terms in the fits which would necessitate non-linear fit- 
ting procedures which usually tend to be quite unstable. Of course, such a 
conclusion depends on the accuracy of the data and only shows that the cor- 
rection terms are so small that they cannot be resolved within our accuracy. 
This is quite remarkable because in the pure 3D Ising model study of Ref. 67 
for lattices of size L = 8 to 128 (using a different FSS technique) conflu- 
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ent corrections to scaling oc with uo = 0.87(9) could be resolved. One 
possible explanation is that random lattices with large average coordination 
number and preserved rotational invariance indeed reach the infinite-volume 
limit earlier than their regular counterparts. This is quite conceivable since 
the magnitude of the correction term does not only depend on the exponent 
but also on its amplitude, and the latter is a non-universal quantity which 
does definitely depend on the lattice structure. 

Having now 28 estimates of 1/u from fits with Q > 0.15 we proceeded as 
follows. We computed arithmetic and error weighted averages for the four 
subgroups of estimates and finally over the total set of estimates. The main 
results are collected in Table 5. Because of the neglected cross-correlations 
between the fit results, in particular the error estimate of the weighted av- 
erage should be taken with some care. As our best estimates we therefore 
quote throughout this paper in all tables in the lines labeled "final" always 
the weighted mean and quite conservatively the smallest error among all 
available fits (making thus the plausible assumption that the error of the 
weighted average is never larger than the error of the most accurate fit re- 
sult that contributes to this average). As a general trend we see in Table 5 
that the fits of d ln[(|m| p )] av /<iK" yield more accurate results than those of 
dU2 P /dK. We also observe, however, that the partial averages over the fits 
of dU2 P /dK,p = 1, 2 and rfln[(|m| p )] av /<iK', p = 1, 2, are only marginally con- 
sistent, even though in all cases the goodness-of-fit was high. We thus have 
no reason based on statistical arguments to favor one or the other group of 
fits and therefore take as our final value the weighted average over all 28 
estimates which yields 

l/v = 1.5875 ±0.0012, v = 0.62992 ± 0.00048, (21) 

with the minimal error coming from the fit of d \n[{\m\ 2 )} av /dK at the max- 
imum locations of d[(\m\)] av /dK. 

If we only average the results of the fits of the maxima of d ln[(]m|)] av / dK, 
dln[(m 2 )] aY /dK, dU 2 /dK, and dU^/dK, we obtain basically the same final 
estimate with a slightly larger error bar: 

l/v = 1.5878 ±0.0015, v = 0.62980 ± 0.00060. (22) 

To give also a visual impression of the quality of these fits, they are shown 
in Fig. 3. This reconfirms the weighted average (21) over all 28 fits, and our 
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Table 5: Fit results for the critical exponent \ jv. 
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final estimate for the correlation length exponent is thus 

v = 0.6299 ± 0.0005. (23) 

For comparison, for Ising models on regular simple cubic (SC) lattices Fer- 
renberg and Landau 68 obtained v = 0.6289(8), Blote et al. m concluded that 
v = 0.6301(8), Ballesteros et al. 67 found 70 v = 0.6294(5) [5], and in a recent 
study of the 4 lattice field theory Hasenbusch 71 estimated v = 0.6296(3) [4]. 
Within error bars all three estimates are in perfect agreement with our result 
(23) for random lattices. 

B. Critical coupling K c : Having determined the critical exponent z/, it 
is straightforward to obtain estimates of the critical coupling K c from linear 
least-squares fits to 

K maxi = K C + a l L~ 1 ^, (24) 

where K max . are the seven pseudo-transition points discussed earlier. Here 
we found a significant improvement of the quality of the fits if the smallest 
lattice size with N = 2 000 was excluded. This can also be inspected visually 
in Fig. 4, where the data and fits are shown. We see a systematic trend that 
the N = 2 000 data lie a little bit too low. In Table 6 we therefore display 
the fit results over the six lattice sizes N = 4 000 — 128 000. By using the 
same averaging procedure as before we arrive at the final estimate 

K c = 0.0724 249 ± 0.0000 040. (25) 

Of course, in principle this estimate is biased by our estimate of v. We have 
checked, however, that the dependence on v is extremely weak. If we repeat 
the fits with \jv — 1.5875 ± 3eiA,, where e\i v = 0.0012 is the error on the 
estimate of \jv in (21), we obtain a variation in the estimate for K c by only 
±2 in the last digit, which is much smaller than the statistical error in (25). 

C. Critical exponent 7: The exponent ratio 7/z/ can be obtained from 
fits to the FSS behavior (16) of the susceptibility. By monitoring the quality 
of the fits we decided to discard the N = 2 000 data for the if maxc and 
K maXdm/dK sequences (which led to Q values of 0.05 and 0.02, respectively). 
The fits collected in Table 7 then all have Q > 0.25. The final result is 

7/1/ = 1.9576 ±0.0013, (26) 
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Table 6: Fit results for the critical coupling K c , using the FSS ansatz K maXi = 
K c + aiL~ x l v and our best estimate of \ jv = 1.5875(12). The fit range is 
always from N = 4 000 to 128 000. 
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which should be compared with the estimates for regular SC lattices of 7/V = 
1.970(14) in Ref., 68 7/1/ = 1.9630(30) in Ref., 69 7/// = 1.9626(6)[6] in Ref., 67 
and 7/1/ = 1.9642(4) [5] in Ref.. 71 

For the exponent 77, the estimate (26) implies 

77 = 2 - 7/1/ = 0.0424 ± 0.0013, (27) 

and, by using our value (21) for l/u, we derive 

7 = 1.2332 ±0.0018. (28) 



D. Critical exponent /3: The exponent ratio /3/v can be either obtained 
from the FSS behavior of [(|w|)] av or d^m^]^ / dK , Eqs. (15) or (17). In the 
first case, the sequences K maXx and K ma _ Xdm/dK yield poor Q values (< 0.01) 
if the N = 2 000 are included in the fits. If we discard the smallest lattice in 
these two cases, all fits shown in Table 7 are characterized by Q > 0.10. The 
final estimate is then 

P/v = 0.51587 ± 0.00082, (29) 
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Table 7: Fit results for the critical exponents 7/1/, (3/v, and (1 — f3)/v. The 
superscripts * and # at the Q values indicate that these fits start at N = 4 000 
and N = 8 000, respectively. The other fits use all data from N = 2 000 to 
128 000. 
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1.9548(38) 
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1.06822(84) 
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and, by using our estimate for \ jv in (21), 

(3 = 0.32498 ± 0.00077. (30) 

If hyperscaling is valid, the estimate (29) would imply ^ /u = D — 2(3/ v = 
1.9683(17), which however turns out to be only barely consistent with the 
direct estimate (26) oi^/v. 

The FSS of d[{\m\)] av / dK is less well behaved. Here we had to discard 
for the K maxc , K mSLXdm/dK , K mSLXd[nm/dK , and K mSLXd[nm2/dK sequences both the 
N = 2 000 and N = 4 000 data in order to guarantee that all fits entering 
the average have a goodness-of-fit parameter Q > 0.15. We then obtain 

(1 -(3)/v= 1.0680 ± 0.0014, (31) 

and by inserting the estimate (21) for 1/v, 

(3/v = 0.5194 ±0.0026, (32) 

and 

(3 = 0.3272 ±0.0014. (33) 

Recent MC estimates for regular SC lattices are (3/v — 0.518(7) in Ref. 68 
and (3/v = 0.5185(15) in Ref.. 69 
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E. Critical exponent a: Due to the regular background term C rcg in 
the FSS behavior (14), the specific heat is usually among the most difficult 
quantities to analyse. 72 We tried non-linear fits to the ansatz C = C rcg + 
aL a / u , but for most sequences of pseudo-transition points the errors on the 
parameters of this fit turned out to be large. We therefore fixed the exponent 
ajv at the value one would expect if hyperscaling is valid, 

a/v = 2/v — D — 0.1750 ± 0.0024, (34) 
a = 2-/^ = 0.1102 ±0.0015, (35) 

and tested if linear two-parameter fits yield acceptable goodness-of-fit values. 
The results are shown in Fig. 5. We see that over the whole range of lattice 
sizes the expected linear behavior is satisfied. The quantitative analysis 
reveals some deviations for the two smallest lattice sizes, but for the fits 
starting with N = 8 000 we obtained for all seven sequences of pseudo- 
transition points goodness-of-fits parameters Q > 0.5. 

F. Binder parameters U2 and U4: It is well known 73 that the U2 P (K) 
curves for different lattice sizes L should intersect around (K c , U 2p ) with 
slopes U 2p = dU 2p /dK oc L 1 ^, where U 2p is the (weakly universal) "renor- 
malized charge". In Fig. 6 we show U 2 and U4 as a function of N = L 3 for 
5 K- values around K c « 0.072 425. At our best estimate of K c , both cu- 
mulants seem indeed to be almost independent of the lattice size. Taking as 
final estimate the weighted mean value (i.e., a least-squares fit to a constant) 
over the results for N = 8 000 - 128 000, we obtain 

U; = 0.58706 ± 0.00044, (36) 
Ul = 0.4647 ±0.0012. (37) 

The variation due to the uncertainty in K c is about twice the statistical 
error at fixed K (0.00080 for U 2 and 0.0020 for U4). For comparison, for 
the standard nearest-neighbor Ising model on a SC lattice, Ferrenberg and 
Landau 68 estimated Ul ~ 0.47, by combining results for three different spin 
models belonging to the Ising universality class, Blote et al. 69 derived C/| = 
0.4652(4), and Ballesteros et al. obtained U* A = 0.4656(4)[4]. For the <p A 
lattice field theory Hasenbusch 71 extracted Ul w 0.46555(9). 
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6 Concluding remarks 



We have performed a detailed finite-size scaling analysis of single-cluster 
Monte Carlo simulations of the Ising model on three-dimensional Poissonian 
random lattices of Voronoi/Delaunay type. At first sight our use of different 
quantities to estimate the same critical exponent might appear redundant, 
since the various estimates are, of course, not independent in a statistical 
sense. Their consistency, however, gives confidence that corrections to the 
asymptotic scaling behavior are very small and can safely be neglected. Our 
estimates for the exponents u, fi/v, and 7/V are all consistent with the best 
numerical estimates for the three-dimensional Ising model and 4 field theory 
on regular lattices - at a very high level of accuracy which is comparable with 
the best estimates coming from field theoretical techniques, cf. Table 8. 74 
While our exponent ratio 7/V would also be compatible with recent estimates 
for the 3D Ising disordered fixed point, our estimate for U\ is more consistent 
with the pure Ising model estimates. The cleanest result yields the critical 
exponent u, where our result agrees within error bars with all previously 
derived estimates for the pure model but is clearly incompatible with the 
disordered fixed point value. We thus obtain strong evidence that, for the 

Table 8: Recent estimates of critical parameters of the pure and disorderd 3D 
Ising model (SC = simple cubic lattice, RG = renormalization group, SD = 
site- dilution, RIM = random Ising model). 
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considered lattice sizes up to iV = 128 000 ~ 50 3 , the Ising model on three- 
dimensional Poissonian random lattices of Voronoi/Delaunay type behaves 
effectively as on regular lattices. 

Of course, we cannot exclude the possibility that on much larger length 
scales (lattice sizes) the scaling behavior may change. Such a late crossover 
is conceivable in the case of weak disorder, where the asymptotic critical 
behavior governed by a "disordered" fixed point may show up only in the 
extremely close vicinity of criticality, that is at extremely large system sizes in 
a finite-size scaling analysis. Even though the qualitative scaling behavior is 
expected to be universal, quantitative properties of the crossover point such 
as its location should depend on the strength of the disorder via non-universal 
amplitudes. In order to obtain for the random lattices a rough estimate 
of the strength S of the local connectivity disorder we have computed the 
relative variance of the local coordination numbers which may be viewed as 
a measure for the size of effective temperature variations over the lattice. 
From the probability density P(q) displayed in Fig. 2 it is straightforward to 
obtain 

S = (q - q) 2 /f = 0.0461, (38) 

withg = 2+487r 2 /35 = 15.5354 Similarly, for two-dimensional Poissonian 

Voronoi/Delaunay random lattices one finds S = 0.0491 with q — 6. The 
relative variance (38) can be compared with the fluctuations of the number 
B of active bonds per site in bond-diluted models. Here B follows a binomial 
distribution, P(B) = (^ 2 ^p B (l — p) 2D ~ B , where D is the dimension and p 
denotes the probability for a bond to be active (such that p — 1 corresponds 
to the pure model), and one obtains 

1 1 -p 



S = (B - BY/B = (39) 

v ' 1 2D p 1 K ' 

with B = 2Dp. By equating Eqs. (38) and (39) and solving for the dilution 
parameter p one can thus determine an associated bond-dilution model with 
the same local disorder fluctuations as for the random lattices. For the three- 
dimensional case this yields p = 0.7834, and in two dimensions one finds 
p = 0.8358. In the terminology of three-dimensional bond-diluted Ising 7 " 9 
and g-state Potts 17 ' 19 models such a value of p belongs to the weak dilution 
regime where some influence of the disordered fixed point can be observed, 
but it is still difficult to clearly disentangle it. For site-diluted models the 
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corresponding p-value is presumably higher, in particular for weak dilution, 
since all bonds around a vacant site are non-active. In the latter models, of 
course, the dilution parameter p can easily be tuned to study more accessible 
regions. 

In view of the very high quality of our fits based on the leading FSS ansatz 
only we must conclude that in the case of Voronoi/Delaunay random lattices 
very much larger system sizes would be necessary to observe the expected 
crossover to the critical behavior associated with the disordered fixed point. 
This was clearly outside the scope of the present study and its computer 
budget which was equivalent to several years of fast workstation CPU time. 
Instead of further increasing the system size, an alternative and more promis- 
ing route for future studies could be a systematic variation of the random 
lattice construction by modifying the Poissonian nature of the site distribu- 
tion such as to achieve larger values of S (corresponding to smaller values 
of p) , or an investigation of the present random lattices coupled to a model 
with a larger critical exponent a where the expected crossover should set in 
for reasonable lattice sizes already for a moderate degree of local disorder. 
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Appendix: Random lattice construction 



The employed algorithm for the random lattice construction works as fol- 
lows. Adapting the method described in Ref., 60 we first draw randomly iV 
sites uniformly distributed in a unit volume, thereby approximating a Pois- 
sonian distribution. For alternative distributions discussed in the literature 
see, e.g., Refs.. 41 ' 59 In the second step we link the sites according to the 
Voronoi/Delaunay prescription. We start by picking the first site that we 
drew and locate all its nearest neighbors, keeping them stored in an array. 
Then we proceed to the second site and search for all its nearest neighbors 
too; once finished with the second site we keep repeating the procedure until 
we have done it for all the sites. Of course, with this method we locate a 
given link twice, but the simplicity of its implementation pays off. 

Starting from a given site, x±, the linking procedure works as follows. Its 
nearest neighbor, X2, is located from within the few hundred sites forming, 
or belonging to, the "cloud-of-neighbors" around X\. We will comment later 
the issue of how to determine a cloud-of-neighbors for a given site. Notice 
that some care must be exercised when approaching the boundaries of the 
lattice to ensure the periodic boundary conditions. Afterwards a third site, 
X3, is searched for, and a triangle is constructed. In order to locate this third 
site, we draw circumferences going through xi, £2, and the few hundred 
sites belonging to the cloud-of-neighbors of x±. We pick as x 3 the site for 
which the radius of the circumference is the smallest. From the triangle 
A(xi, X2, X3) we proceed to locate a fourth site, X4, linked to x\ and build 
a tetrahedron, t(xi, X2, £3, X4). When a triangle still has not been used to 
build a tetrahedron from it, it is termed "active" and a logical flag is turned 
"on" . When already used, it is renamed to "inactive" and the flag is turned 
"off". 

To construct a tetrahedron from a triangle we split the volume in two half 
spaces: one "above" and the other "below" the plane lying on the triangle 
A(xi, #2, xs). Let us suppose that we search for X4 in the half space "above" 
the triangle. In order to determine £4, we draw spheres going through x±, 
X2, X3 and the sites belonging to the cloud-of-neighbors of x\ placed in the 
half space in which we are working. If we happen to find several trial sites 
for which their distance to the circumcenter of the triangle A(xi, X2, £3) is 
smaller than the radius of the circumscribed circle of A(xi, X2, X3), then 
we pick as rr 4 the site for which the radius of the circumscribed sphere of 
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t(xi, x 2 , X3, X4) is the biggest. If, on the contrary, all the trial sites lie at 
a distance from the circumcenter of the triangle A(x±, X2, £3) greater than 
the radius of the circumscribed circle of A(xi,x 2 , x 3 ), then we pick as X4 the 
site for which the radius of the circumscribed sphere of t(x±, x 2 , x 3 , X4) is 
the smallest. From the newly built tetrahedron t(xi, x 2 , £3, X4), we can take 
two "active" triangles, A(x±, x 2 , X4) and A(x\, £3, X4) to continue building 
tetrahedra from triangles, and then triangles from tetrahedra. The closer 
neighbors of a given site x\ are all found when there is no "active" triangle 
left connected to the site. 

When describing how to locate xi's nearest neighbor, x 2 , or how to find X3 
afterwards we emphasized that we only search from within the sites forming a 
cloud-of- neighbors around x\. Its meaning is that before starting the linking 
procedure we set up an array for each site in the lattice containing the sites 
forming its cloud. A given site will belong to the cloud of, say, X\ if it lies 
within a sphere centered in x\. The radius of the sphere is chosen such that, 
on the average, the number of sites within the sphere is three times of an a 
priori upper limit to the maximum number of links that a site is likely to 
have in a finite Voronoi/Delaunay random lattice. To implement an efficient 
search of the sites which will belong to the cloud-of-neighbors of a given site, 
we subdivided the unit volume into smaller boxes. The optimal box size 
should be large enough to ensure that nearest neighbors will be located in 
the same box or at least in one of the 26 surrounding boxes, but small enough 
to minimize the time needed for testing all trial sites in a box. 
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Figure 1: The total CPU time spent in constructing a three-dimensional Pois- 
sonian random lattice according to the Voronoi/Delaunay prescription versus 
the number of sites N . The circles show the fraction of time needed to set up 
the " cloud- of-neighbors " . 
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Figure 2: (a) The probability density P(q) of the coordination numbers 
The average is q = 2 + 487r 2 /35 = 15.5354 .... (b) The same data as in \ 
on a logarithmic scale. 




Q- 



10- 6 - 



33 





Figure 4: FSS fits of the pseudo-transition points K max . with \jv — 1.5875 
fixed, yielding a combined estimate of K c = 0.0724 249(40). 
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Figure 5: FSS behavior of the specific heat, assuming a/u = 2/v — D = 
0.1750. 
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Figure 6: FSS behavior of the magnetic cumulants. The central value of 
K is our best estimate (25) for the inverse critical temperature. For the 
neighboring curves the K values vary by about one statistical error bar. 
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